
#include <cusp/precond/detail/device/bspai.h>
#include <cusp/multiply.h>

namespace cusp
{
namespace precond
{

/* Create with same sparsity pattern as A */
template <typename IndexType, typename ValueType, typename MemorySpace>
template <typename MatrixType, typename SparsityPatternType>
bspai<IndexType, ValueType, MemorySpace>::bspai(const MatrixType& A, const SparsityPatternType& A_BSP)
	: linear_operator<ValueType,MemorySpace, IndexType>(A.num_rows, A.num_cols)
{
	make_bspai(A, A_BSP, A_BSP);
}


/* Create with sparsity pattern P */
template <typename IndexType, typename ValueType, typename MemorySpace>
template <typename MatrixType, typename SparsityPatternType>
bspai<IndexType, ValueType, MemorySpace>::bspai(const MatrixType& A, const SparsityPatternType& A_BSP, const SparsityPatternType& M_BSP)
	: linear_operator<ValueType,MemorySpace, IndexType>(A.num_rows, A.num_cols)
{
	make_spai(A, A_BSP, M_BSP);
}

template <typename IndexType, typename ValueType, typename MemorySpace>
template <typename MatrixType, typename SparsityPatternType>
void bspai<IndexType, ValueType,MemorySpace>::make_bspai(const MatrixType& A, const SparsityPatternType& A_BSP, const SparsityPatternType& M_BSP)
{
    cusp::precond::detail::device::spai_core(A, A_BSP, M_BSP, matrix);
}


template <typename IndexType, typename ValueType, typename MemorySpace>
template <typename VectorType1, typename VectorType2>
void bspai<IndexType, ValueType,MemorySpace>::operator()(const VectorType1& x, VectorType2& y) const
{
	cusp::multiply(matrix, x, y);
}


} // end namespace precond
} // end namespace cusp

